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Abstract 

In  an  effort  to  unify  the  various  phase-field  models  that  have  been  used  to  study 
solidification,  we  have  developed  a class  of  phase-field  models  for  crystallization  of  a 
pure  substance  from  its  melt.  These  models  are  based  on  an  entropy  functional,  as  in 
the  treatment  of  Penrose  and  Fife,  and  are  therefore  thermodynamically  consistent 
inasmuch  as  they  guarantee  spatially  local  positive  entropy  production.  General 
conditions  are  developed  to  ensure  that  the  phase  field  takes  on  constant  values 
in  the  bulk  phases.  Specific  forms  of  a phase-field  function  are  chosen  to  produce 
two  models  that  bear  strong  resemblances  to  the  models  proposed  by  Langer  and 
Kobayashi.  Our  models  contain  additional  nonlinear  functions  of  the  phase  field 
that  are  necessary  to  guarantee  thermodynamic  consistency. 

PACS  numbers:  64.70.Dv,  68. 10. Gw,  81.30.Bx,  82.65.Dp 

Keywords:  non-equilibrium  thermodynamics,  phase  field,  solidification,  free  boundary, 
entropy  functional.  Model  C 


* Consultant,  National  Institute  of  Stcindards  and  Technology 

^Technology  Administration,  U.S.  Department  of  Commerce,  Washington  D.C. 


1 


1.  Introduction 


The  classical  approach  to  the  modeling  of  first  order  phase  transformations  involves  track- 
ing of  the  free  boundary  that  separates  the  growing  phase  from  the  parent  phase.  This 
requires  the  solution  of  a formidable  free  boundary  problem.  For  example,  for  crystal- 
lization of  a solid  (crystal)  from  a pure  liquid  (melt),  one  must  solve  the  equations  for 
the  temperature  fields  in  both  solid  and  liquid  subject  to  boundary  conditions  for  the 
temperature  and  its  derivatives  on  the  moving  sohd-liquid  interface,  the  free  boundary. 

The  phase-field  model  provides  an  alternative  approach,  according  to  which  a new 
variable,  the  phase  field  (f),  is  introduced  to  keep  track  of  the  phase,  taking  on  constant 
values  indicative  of  each  of  the  bulk  phases  and  making  a transition  between  these  values 
over  a thin  transition  layer  that  plays  the  role  of  the  classically  sharp  interface.  The  phase 
field  (f)  is  governed  by  a partial  differential  equation  that  guarantees  (in  the  asymptotic 
limit  of  a suitably  thin  transition  layer)  that  the  appropriate  boundary  conditions  at  the 
crystal-melt  interface  are  satisfied.  Moreover,  the  transport  equations  are  modified,  by 
addition  of  terms  that  depend  on  the  phase  field,  to  apply  in  all  of  space  and  to  embody  the 
conservation  conditions  at  the  interface.  One  then  proceeds  to  solve  the  coupled  equations 
for  the  phase  field  and  transport-related  fields  (e  g.,  temperature  and/or  composition) 
without  the  necessity  of  exphcit  tracking  of  the  free  boundary. 

Langer  [1,2]  adapted  Model  C of  Halpern,  Hohenberg,  and  Ma  [3]  to  produce  a phase- 
field  model  that  accounts  for  the  Gibbs-Thomson  equation  at  the  sohd-liquid  interface  for 
sohdification  problems.  Moreover,  a similar  phase-field  model  was  introduced  by  Collins 
and  Levine  [4].  Analytical  properties  of  the  phase- field  equations  resulting  from  hanger’s 
model  have  been  analyzed  extensively  by  Caginalp  et  al.  [5-13],  Gurtin  [14]  and  Fife  and 
GiU  [15, 16].  In  particular,  Caginalp  [11]  has  shown  that  certain  distinguished  limits  of 
the  phase-field  equations  lead  to  various  free  boundary  problems,  including  the  classical 
models  of  sohdification. 

Numerical  methods  for  computation  using  a phase-field  model  were  introduced  by 
Fix  [17]  and  early  numerical  computations  in  one  dimension  were  performed  by  Lin  [18]. 
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Subsequently,  there  have  been  a number  of  computations  in  one  dimension,  including  those 
by  Schofield  and  Oxtoby  [19]  for  planar  geometries  and  by  Caginalp  and  Socolovsky  [12] 
for  planar  and  spherically  symmetric  geometries.  Success  with  numerical  computations  for 
more  complex  shapes  has  been  limited  because  of  the  considerable  computational  power 
needed  to  resolve  the  diverse  length  and  time  scales  intrinsic  to  the  phase-field  model. 
Recently,  Kobayashi  [20]  introduced  a phase-field  model  that  is  particularly  suitable  for 
numerical  computations  of  the  dendritic  solidification  of  a pure  material  from  a supercooled 
melt,  in  two  [21,22]  and  three  [23]  dimensions.  Using  the  phase-field  model  described  here, 
numerical  computations  in  two  dimensions  have  been  carried  out  for  the  dendritic  growth 
of  nickel  by  Wheeler  et  al.  [24]. 

A difficulty  with  phase- field  models,  as  pointed  out  by  Penrose  and  Fife  [25],  lies  in 
their  derivation.  In  particular,  the  governing  equation  for  the  phase  field  is  usually  derived 
from  a free  energy  functional  that  is  apphcable  only  to  an  isothermal  situation;  then  the 
classical  equation  for  the  temperature  field  is  altered  in  an  ad  hoc  manner  to  account  for 
the  hberation  of  latent  heat  by  addition  of  a term  proportional  to  the  time  derivative  of  the 
phase  field.  Penrose  and  Fife,  however,  provided  a framework  from  which  the  equations 
for  the  phase  field  and  the  temperature  can  be  derived  in  a thermodynamically  consistent 
manner  from  a single  entropy  functional. 

In  the  present  paper,  we  employ  the  method  of  Penrose  and  Fife  to  derive  a class  of 
phase-field  models  for  the  crystallization  of  a pure  material  from  its  melt  that  have  the 
following  properties: 

1.  They  are  based  on  an  entropy  functional  and  guarantee  local  positive  entropy  pro- 
duction, 

2.  The  functions  and  parameters  in  the  model  can  be  chosen  to  agree  with  empirical 
data,  such  as  specific  heats  and  surface  tension,  asymptotically  in  the  classical  limit, 

3.  They  satisfy  criteria  to  insure  that  the  phase  field  will  take  on  fixed  values,  which 
we  take  to  be  0 for  the  soHd  and  1 for  the  liquid,  for  the  bulk  phases. 
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4.  They  allow  for  specific  choices  of  the  phase-field  function  that  lead  to  models,  which 
we  call  Model  I and  Model  II,  that  bear  a strong  resemblance  to  the  models  of 
Langer  and  Kobayashi,  respectively,  but  that  display  specific  nonhnearities  necessary 
to  satisfy  properties  (T3).  Thus,  we  provide  a common  basis  for  their  models. 


2.  Thermodynamic  Formulation 


For  purposes  of  deriving  the  partial  differential  equations  for  the  phase-field  model,  we 
consider  a closed  system  of  volume  V in  which  a pure  material  undergoes  a first-order 
phase  transition  between  sohd  and  liquid.  We  introduce  an  order  parameter  </»(x,t),  called 
the  phase  field,  that  indicates  the  phase  of  the  material  at  (x,  i).  We  assume  that  (f)  equals 
zero  in  the  sohd  phase  and  unity  in  the  hquid  phase.  Thin  regions,  in  which  (j)  takes 
values  between  zero  and  unity,  will  be  shown  to  correspond  to  crystal-melt  interfaces. 
For  simphcity,  we  assume  uniform  density  throughout  the  system  and  that  there  is  no 
convection  in  the  hquid;  this  ensures  mechanical  equilibrium. 

We  assume  that  the  internal  energy  of  any  subvolume  V of  F is  represented  by  the 
integral 

^ — j (1) 

According  to  the  first  law  of  thermodynamics,  the  time  rate  of  change  of  S obeys 

8 f q - nda  = 0,  (2) 

J A 

where  A is  the  surface  of  V with  outward  normal  n,  and  q is  flux  of  internal  energy. 
Converting  the  surface  integral  in  Eq.  (2)  to  a volume  integral  by  Gauss’  theorem  and 
equating  the  integrands,  which  is  permissible  since  V is  arbitrary,  leads  to  the  familiar 
equation 

e-|-V-q  = 0.  (3) 


The  entropy  of  any  subvolume  V of  F is  assumed  to  be  represented  by  the  functional 


dv. 


(4) 
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where  5(e,  is  an  entropy  density  (an  extension  of  classical  thermodynamics)  that  would 
apply  to  a uniform  system  of  internal  energy  density  e and  phase  field  and  e is  a 
constant.  The  second  term  in  the  integrand  is  a gradient  entropy  term  analogous  to  the 
Landau-Ginzburg  or  Cahn-Hilliard  gradient  energy  term  in  the  free  energy,  cf.  Eq.  (37). 
Differentiation  of  Eq.  (4)  with  respect  to  time  gives 


where  we  have  substituted  for  e from  Eq.  (3)  and  used  Gauss’  theorem  to  integrate  by 
parts. 

The  second  law  of  thermodynamics  requires  the  entropy  production  to  be  positive  in 
any  subvolume  V of  V.  This  entropy  production  can  be  calculated  by  subtracting  from  S 
the  entropy  flux  through  A,  which  leads  to  the  relation 


5 + J I y • n|  da  > 0, 


(6) 


where  T is  the  absolute  temperature.  In  Eq.  (6),  q/T  is  the  familiar  entropy  flux  due  to 
heat  flow;  the  quantity  is  an  entropy  flux  related  to  changes  in  the  phase  field  at 

the  boundary  of  the  sub  volume  V.  The  volumetric  entropy  production  and  the  entropy 
flux  associated  with  0 will  be  related  to  their  counterparts  for  a sharp-interface  model  in 
the  discussion  section. 

Substitution  of  Eq.  (5b)  into  Eq.  (6)  then  leads  to 


(a 


qv  - + 


-f 


^ > dv  > 0 


(7) 


where  we  have  used 


de  = Tds  + f — ) d0 
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de 

d(f) 


(8) 


to  identify  {dsjde)^  = (1/T).  We  can  therefore  guarantee  local  positive  entropy  produc- 
tion by  choosing 


q = MtV  ^ , 


(9) 

(10) 


where  Mt  and  r are  positive.  This  choice  allows  the  integrand  in  Eq.  (7)  to  be  reduced 
to  a sum  of  squares:  Iq^/Mr  + The  governing  equations  therefore  become 

1 


e = - V- 

for  the  internal  energy  density,  and 


TJl 


(11) 


(12) 


for  the  phase  field. 

In  order  to  identify  the  derivative  in  Eq.  (12),  it  is  convenient  to  work  with  the 
Helmholtz  free  energy  density 

f = e-Ts, 


whose  differential,  when  combined  with  Eq.  (8),  yields 


and 


ld[SIT]\ 

V dT  A 


The  latter  can  be  integrated  at  constant  (f)  to  give 


f(T,4>)  = T 


+ G(4,) 


(13) 


(14) 


(15) 


(16) 


where  G{(f))  is  an  undetermined  function  of  ^ and  Tm  is  the  melting  point.  To  proceed 
further,  we  assume  that  the  internal  energy  density  has  the  form 


e = es{T)  + p{ct>)L{T)  = e^iT)  -f  - l]T(r), 


(17) 
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where  es{T)  and  eL{T)  are  the  classical  internal  energy  densities  of  the  solid  and  liquid 
phases,  respectively,  and  L{T)  = gl^T)  — es{T)  is  their  difference.  Here  p{(f))  satisfies 
p(0)  = 0 and  p(l)  ==  1.  We  note  that  Lq  = L{Tm)  is  the  latent  heat  of  fusion;  whereas 
the  variation  of  L{T)  > 0 with  temperature  allows  the  specific  heats  of  solid  and  liquid 
to  depend  differently  on  temperature,  as  illustrated  in  Figure  1.  From  Eq.  (16),  the  free 
energy  density  corresponding  to  this  internal  energy  density  is 


f(TA)  = T\-j^ 


Tm  C 


where 


Q[T)  = / 


- M*)  - ilQ(r)  + GW 

T L{0 


■dC. 


(18) 


(19) 


ITu  0 

We  note  that  Q(T)  is  monotonic  increasing  with  respect  to  its  argument,  T,  and  Q{Tm)  — 
0.  From  Eq.  (14)  and  Eq.  (18)  we  have 


= -TQiT)p'W  + TG'W- 


(20) 


Thus,  the  governing  equations,  Eq.  (11)  and  Eq.  (12),  become 


eL{T)+mL{T)  + \p{<j>)  - l]L{T)  = -V 


Afy  V 


(?)]> 


(21) 


and 


= QiT)p’W  - G'W  + 


(22) 


We  have  assumed  that  ^ = 0 and  (j)  = I correspond  to  bulk  sofid  and  liquid  states  for 
all  values  of  the  temperature.  This  consideration  requires  / to  have,  for  all  values  of  T, 
local  minima  with  respect  to  (/>  ai  ^ = 0 and  (f)  = 1.  We  also  require  / to  be  continuous 
with  respect  to  the  temperature  at  the  melting  point;  therefore 


/(Tm,0) 

= /(^m,!), 

(23) 

d(l> 

<^=0,1 

- 0, 

(24) 

av 

<f>=0,l 

> 0. 

(25) 

These  conditions  require  that 


G(0)  = G(l). 


(26) 

(27) 

(28) 


[G'W-v'iWinU^o,.  = 0. 

[G"(.^)-p"(^)Q(T)l|^^„^  > 0. 

Accordingly  we  choose  G{(f))  to  be  a symmetric  ‘double  well’  potential  with  minima  at 
zero  and  unity  of  the  form 

G((»  = (29) 

where 

9{4>)  = <^^(1  - <?^)^  (30) 

and  a is  a positive  constant;  this  choice  implies  the  convention  /{Tm,  0)  = /{Tm,  1)  = 0.  In 
the  next  section  we  explore  two  possible  choices  for  p(^)  which  lead  to  phase- field  models 
that  may  be  compared  to  those  of  Langer  and  Kobayashi. 

We  continue  with  the  general  development  of  the  phase-field  equations  by  further 
specifying  the  thermodynamic  description  of  the  system.  For  a pure  material,  one  is 
mostly  interested  in  studying  morphological  instabilities  and  subsequent  dendritic  growth 
that  can  occur  in  sohdification  from  an  undercooled  melt,  where  heat  conduction  takes 
place  mostly  in  the  liquid.  Therefore,  in  order  to  obtain  a tractable  problem,  we  assume 
the  energy  density  in  the  liquid  to  be  a linear  function  of  temperature  of  the  form 

eL{T)  = eL{TM)  + c{T-TM),  (31) 

where  c is  a constant  specific  heat,  which  will  lead  to  the  classical  heat  equation  in  the 
bulk  liquid.  Eq.  (21)  and  Eq.  (22)  therefore  become 

{c  + (p(^)  - 1)L'{T)}  ^ + L{TW)^  = W^T,  (32) 

and 

= Q(T)p'(«  - (33) 

where  we  have  set  Mt  = kT^.  The  thermal  conductivity  k is  assumed  to  be  constant. 
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At  this  point  in  the  development  of  the  equations,  the  parameters  a and  e can  be 
related  to  the  interface  thickness  and  surface  free  energy  by  examining  a one- dimensional 
solution  of  Eq.  (33)  under  equihbrium  conditions,  T = Tm,  for  which  it  becomes 


2(P(t>  1 //  /\ 

£ — - 7-S  (^)  = 0, 


dx^  4a‘ 

with  boundary  conditions  >0asx— > — oo,  and  > 1 as  x 

found  to  be 


tanh 


X 


+ 1 


(34) 

-}-oo.  The  solution  is 

(35) 


where  we  have  chosen  a constant  of  integration  to  locate  the  interface  at  x ==  0.  From  this 
solution  it  is  clear  that  the  interface  has  characteristic  thickness 


S = \/ai 


(36) 


To  calculate  the  surface  free  energy  associated  with  the  one- dimensional  solution,  we 
employ  the  Helmholtz  free  energy  functional  at  temperature  Tm,  which  from  Eq.  (4)  and 
Eq.  (1)  is  given  by 


= S-TmS=  f 
Jv 


p2'T> 


dv. 


(37) 


The  corresponding  surface  free  energy  per  unit  area,  which  is  a surface  excess  quantity 
because  /(Tm,  ^)  = 0 in  the  bulk,  is 


J — oo  Y dx  J 


This  may  be  evaluated  to  yield 


y/26Ti 


M 


12a 


(38) 


(39) 


Hence,  the  constant  a is  related  to  the  material  parameters  and  the  interface  thickness  6. 

In  order  to  proceed  further,  it  is  advantageous  to  nondimensionahze  the  governing 
equations.  We  introduce  a length  scale  w representative  of  the  size  of  the  domain,  and 
use  the  thermal  diffusion  time  w^jK,  where  k is  the  assumed  constant  thermal  diffusivity 
of  the  liquid,  as  the  reference  time  scale.  Employing  these  scales,  we  write 


X = 


t = 


t 


w 


w'^f  k' 


and  u = 


T-Tm 

Tm-To’ 
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where  To  is  an  appropriate  reference  temperature.  Here,  a tilde  is  used  to  denote  the 
dimensionless  equivalent  of  a dimensional  quantity.  The  governing  equations  become 

L'{u)' 


1 + \p(ii>)  - 1]- 


du  L{u)  ,,  -2 


+ eaSQ{u)p'{(j))  - 

m ot  4 


where 


Q{'^)  = I T- 

Jo  h 


L{u)  = L/Lo, 

m 


and 


a = 


v^u;[To]^ 

12caTM 


S = 


+ ( 

c{Tm  - To) 


A, 


Lo 


. 6 
^ j 

w 


m = 


6V28a 

tkTm 


(40) 

(41) 

(42) 

(43) 

(44) 


Below,  we  determine  m explicitly  in  terms  of  the  interface  kinetic  coefficient  and  known 
material  parameters.  We  shall  proceed  to  discuss  Eqns.  (40-44)  for  two  different  choices 
of  p(0)  given  below. 


3.  Phase-Field  Models  I and  II 

In  order  to  proceed  in  the  development  of  the  phase- field  model,  the  function  p(<^)  which 
appears  in  the  energy  density  must  be  specified.  We  consider  two  different  polynomial 
forms  for  the  function  p{(f>)  which  result  in  different  phase-field  models. 

Model  I:  We  choose 

= I Vw/  = (^5) 

Jo  ^(C)“C 

which  satisfies  the  required  normalization  that  p(0)  = 0 and  p(l)  = 1 as  well  as  the 
conditions  Eq.  (27)  and  Eq.  (28),  independent  of  Q{T),  because  = p"{<t>)  = 0 at 

^ = 0 and  1.  This  choice  will  lead  to  a phase-field  model  that  we  will  compare  with  the 
one  proposed  by  Langer  [2]. 

Model  II:  We  choose 

p{4,)  = 4'\i  - 24),  (46) 
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which  satisfies  the  conditions  that  p(0)  = 0 and  p(l)  = 1,  and  p'(O)  = p'(l)  = 0;  however, 
because  G"(0)  - p"(0)Q(r)  = l/(2a)  - 6Q(T)  and  G"{1)  - p"(l)Q(T)  = l/(2a)  + 6g(T), 
we  must  require  the  additional  condition  that  \\2aQ{T)\  < 1 in  order  to  satisfy  Eq.  (28). 
This  choice,  as  shown  below,  leads  to  a model  similar  to  that  proposed  by  Kobayashi  [20]. 


3.1.  Model  I 


Eq.  (45)  apphes  and  there  are  no  restrictions  on  Q{T).  Therefore,  for  simphcity,  we  choose 
the  specific  heat  of  the  sohd  to  equal  the  constant  specific  heat,  c,  of  the  liquid,  in  which 
case  L{T)  = To,  which  results  in  L{u)  = 1.  We  also  assume  that  \Tm  — To|  <<  Tm  and  so 
we  make  the  approximation  Q{u)  = u.  Thus,  Eq.  (40)  and  Eq.  (41)  become 

du  30g{(l))  d<p 


dl^ 


S di 


= vx 


— ^ + ^0g{(j))eaSu  - jg'{(f>). 

m ot  4 


(47) 


(48) 


The  dimensionless  parameters  S and  a can  be  specified  in  terms  of  ly,  {Tm  — To),  and  the 
material  parameters.  However,  the  parameters  e and  m are  still  undetermined  because 
they  are  related  to  the  unknown  quantities  S and  r.  We  expect  the  parameter  8 to  be 
much  less  than  w.  Hence,  we  fix  the  parameter  m by  conducting  an  asymptotic  analysis  of 
the  phase-field  equations  in  the  limit  e — > 0 with  5,  a,  and  m of  order  one.  In  this  hmit,  a 
free  boundary  problem  is  recovered.  The  details  are  analogous  to  those  given  by  Caginalp 
[11]  and  so  for  brevity  we  just  give  the  results.  The  free  boundary  problem  satisfied  by 
the  leading  order  temperature  is 

| = VV  (49) 

with  interfacial  boundary  conditions 


u 


and 


du 


dn 


1 Liquid 
Sohd 


■5’ 


(50) 


(51) 
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where  f = aTml[yjLo{TM  — ^o)]  is  the  dimensionless  capillary  length,  vj  is  the  dimen- 
sionless normal  velocity  of  the  interface  into  the  liquid  and  /C  is  the  dimensionless  interfa- 
cial curvature  (measured  positive  for  convex  projections  into  the  liquid).  The  interfacial 
boundary  condition  Eq.  (50)  reduces  to  the  Gibbs-Thomson  equation  when  vj  = 0.  The 
additional  term  incorporates  the  effect  of  interface  kinetics.  The  dimensionless  interface 
kinetic  coefficient  is  therefore  represented  by  m;  it  may  be  related  to  the  corresponding 
dimensional  quantity  by  comparison  of  Eq.  (50)  with  its  usual  dimensional  form; 

T = T„-  - -V,,  (52) 

Lq  fj- 

which  gives 

fiaTM 
m = — — 
kLq 

In  view  of  Eq.  (44),  this  means  that  r is  proportional  to  6,  which  makes  m independent 
of  e as  implied  by  Eq.  (41).  The  only  remaining  undetermined  parameter  is  therefore  e. 
In  practice,  e should  be  chosen  small  enough  to  validate  the  asymptotics  but  large  enough 
to  make  computations  practical.  This  completes  the  determination  of  the  dimensionless 
phase- field  parameters  in  terms  of  the  physical  parameters. 

The  advantages  of  the  phase-field  model  given  by  Eq.  (47)  and  Eq.  (48)  are  two-fold. 
First,  the  model  is  consistent  with  the  ideas  of  irreversible  thermodynamics  because  it  was 
derived  from  an  entropy  functional,  which  is  the  appropriate  thermodynamic  potential  in 
this  nonisothermal  situation.  Second,  it  is  constructed  so  that  the  states  ^ = 0 and  (f)  = 1 
correspond  to  the  bulk  sohd  and  hquid  phases,  respectively,  independent  of  the  temper- 
ature. This  has  the  advantage  that  latent  heat  of  the  correct  total  amount  is  hberated 
only  in  the  interfacial  region.  The  model  proposed  by  Langer  [1,2]  can  be  obtained  by 
replacing  the  term  30p(^),  which  appears  in  both  Eq.  (47)  and  Eq.  (48),  by  unity,  which 
is  its  average  value  over  the  range  0 < ^ < 1,  while  leaving  the  term  in  Eq.  (48) 
as  calculated  from  Eq.  (30).  This  “hneaxization”  would  violate  Eq.  (27)  and  Eq.  (28). 
Thus,  in  the  Langer  model,  the  bulk  solid  and  liquid  phases  do  not  correspond  to  constant 
values  of  (f)  independent  of  the  temperature.  In  a recent  model  developed  by  Caginalp  and 
Chen  [13],  a function  /(0)  was  introduced  in  the  free  energy  functional  (their  equation  1.6) 
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to  identify  the  bulk  phases  precisely  by  (f)  = ±1,  which  is  equivalent  to  maintaining  the 
term  30^(^)  in  Eq.  (48),  but  their  model  does  not  have  the  thermodynamically-consistent 
30g{(l))  term  in  Eq.  (47). 


3.2.  Model  II 


Eq.  (46)  apphes,  but  there  is  the  restriction  \12aQ{T)\  < 1.  Substituting  Eq.  (46)  for  p(^) 
and  Eq.  (30)  for  into  Eq.  (40)  and  Eq.  (41)  yields  the  following  set  of  phase-field 
equations: 


d(t> 


dL[u)  6L{u) 


di 


u. 


(54) 


— ^ + il>{i  - <l>)  (4  - ^ + SeaSQ{u)^  . (55) 

Recently,  Kobayashi  [20-23]  proposed  a phase-field  model  based  on  the  following  pair 
of  dimensionless  equations: 

fill  1 r)rh  ~ . 

(56) 


du  1 dcj)  - , 


~ (f}{l  — (f))  9'rctan(i/u) 


(57) 


where  /3  and  u are  positive  constants  with  /3  < 1.  The  phase- field  equation  in  Model  II, 
Eq.  (55),  is  identical  to  that  suggested  by  Kobayashi  provided  that 

/? 


Q{u)  = 


direaS 


arctan(i/u). 


(58) 


Note  that  our  restriction  |12aQ(r)|  < 1 is  equivalent  to  Kobayashi’s  condition  /?  < 1.  The 
derivative  of  the  function  Q(u)  is  related  to  the  function  L(u)  by  Eq.  (43),  which  gives 


l(u)  = 


[i  + (^) 


u 


(59) 


1 4-  (i/u)2 

and  ^ = 67rea5/i/.  The  function  L{u)  will  have  a strong  dependence  on  u unless  v is 
chosen  to  be  very  small.  This  dependence  on  u requires  the  specific  heat,  C5,  of  the  solid 
to  have  a definite  dependence  on  temperature,  which,  from  the  derivative  of  the  expression 
es{T)  = ^l{T)  — L(T),  Eq.  (31),  and  Eq.  (42),  is  represented  by 


cs  = c 


1 - |i'(>') 


(60) 
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The  quantity  —L'{u)  is  plotted  in  Figure  2 for  the  case  u = 10  (used  by  Kobayashi)  and 
(Tm  — To)ITm  = 0.2.  We  note  that  the  corresponding  specific  heat  will  vary  considerably, 
even  near  the  melting  point,  u = 0,  and  can  possibly  become  negative.  For  sohdification 
from  a supercooled  melt,  the  temperature  in  the  solid  will  typically  be  close  to  the  melting 
point,  with  temperature  gradients  much  less  than  those  in  the  fiquid.  Nevertheless,  there 
might  be  problems  in  using  Kobayashi’s  model  for  such  large  values  of  v. 

Thus,  in  order  to  recover  the  same  governing  equation  for  the  phase  field  as  Kobayashi, 
within  a thermodynamically-consistent  framework,  it  is  necessary  for  the  specific  heat  of 
the  solid  to  have  the  functional  dependence  on  temperature  given  by  Eq.  (60).  However, 
the  self-consistent  equation  for  the  temperature  will  then  be  Eq.  (54)  instead  of  Kobayashi’s 
Eq.  (56).  We  first  compare  the  terms  containing  d(j)ldt\  they  will  only  be  important  in  the 
interfacial  region  where  u ~ 0 so  L{u)  ~ 1.  Thus,  replacing  60(1—0)  by  its  average  value  of 
1 over  the  range  0 < 0 < 1 reproduces  Kobayashi’s  corresponding  term.  This  is  similar  to 
the  “linearization”  required  to  recover  hanger’s  model  from  Model  I.  The  term  in  Eq.  (54) 
containing  dLjdt,  however,  is  completely  absent  in  Kobayashi’s  Eq.  (56).  It  vanishes  in 
the  hquid  where  0=1,  but  not  in  the  solid  where  0 = 0,  and  arises  because  of  the  variable 
specific  heat  necessary  to  make  Kobayashi’s  model  thermodynamically  consistent.  In  fact, 
for  0 ==  0,  Eq.  (60)  can  be  used  to  combine  the  first  two  terms  on  the  left-hand  side  of 
Eq.  (54)  to  yield 


C5  du  -2 
C Ot 


(61) 


which  is  in  complete  agreement  with  the  classical  equation  in  the  solid  phase. 

An  alternative  approach  to  comparing  our  model  to  that  of  Kobayashi  would  be  to 
choose  C5,  and  p(0)  to  achieve  the  same  equation  for  the  temperature.  This  would  require 
cs  = c and  p(0)  ==  0,  but  then  G'(0)  — p'{(f))Q{T)  = (j'(0)  — Lq  [(I/Tm)  — (1/T)],  and 
Eq.  (27)  would  require  G(0)  and  G(l)  to  depend  on  temperature,  which  is  a contradiction. 
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4.  Discussion 


We  have  employed  the  method  of  Penrose  and  Fife,  based  on  an  entropy  functional,  to 
derive  a class  of  phase-field  models  for  the  crystalhzation  of  a pure  material  from  its  melt. 
This  method  provides  a Lyapunov  functional  [25]  which  may  facihtate  the  development 
of  numerical  solution  methods.  The  parameters  of  the  model  are  selected  by  passing  to 
the  asymptotic  limit  of  a sharp  interface  and  comparing  with  the  classical  model.  Criteria 
are  developed  to  ensure  that  the  phase  field  will  take  on  fixed  values,  independent  of 
temperature,  for  bulk  solid  and  hquid.  Specific  choices  of  the  phase-field  function  are  made 
to  yield  models,  which  we  call  Model  I and  Model  II,  that  bear  strong  resemblances  to  the 
models  of  Langer  and  Kobayashi,  respectively.  However,  Models  I and  II  contain  nonlinear 
functions  of  the  phase  field  that  are  absent  in  the  models  of  Langer  and  Kobayashi,  but 
are  necessary  to  guarantee  thermodynamic  consistency.  Model  I differs  from  Langer’s 
model  only  in  the  interfacial  region.  Model  II  differs  from  Kobayashi’s  model  both  in  the 
interfacial  region  and  in  the  bulk  sohd;  in  particular,  the  equation  for  the  temperature 
that  is  thermodynamically  consistent  with  Kobayashi’s  equation  for  the  phase  field  has  a 
specific  heat  in  the  solid  phase  that  is  a rather  strong  function  of  temperature,  even  for 
temperatures  close  to  the  melting  point. 

These  models  also  guarantee  positive  entropy  production  locally  va.  space.  The  volumet- 
ric entropy  production  due  to  heat  flow,  IqP/Mt,  and  the  entropy  flux,  q/T,  are  classical 
and  require  no  further  discussion.  The  corresponding  quantities,  t(^)^  and  require 

further  clarification,  which  we  proceed  to  give  by  relating  them  to  analogous  quantities  for 
a sharp-interface  model.  The  volumetric  entropy  production  will  contribute  only  in 

volume  elements  dv  in  which  (f>  differs  appreciably  from  zero.  These  elements  will  corre- 
spond to  regions  through  which  a solid-hquid  interface  is  passing,  and  for  which  (j)  rKj  vilv 
and  dv  ~ 7/da/,  where  daj  is  an  element  of  area  on  the  interface.  Here  we  have  taken 
T)  = 6\/26  as  a more  precise  representation  of  the  interface  thickness;  from  Eq.  (35),  rj 
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represents  the  distance  over  which  <f>  varies  from  roughly  0.05  to  0.95.  Thus, 


dv  ~ —vj  daj  = ■■  da/, 


(62) 


V ' f^Tlj 

where  Eqs.  (44)  and  (53)  have  been  used.  The  quantity  [LQ/{fj,T^)]vj  turns  out  to  be 
the  entropy  production  per  unit  area  associated  with  the  movement  of  a sharp  interface 
for  a hnear  kinetic  law  of  the  form  represented  by  Eq.  (52),  as  can  be  shown  using  the 
methodology  of  Gurtin  [26]. 

Similarly,  the  entropy  flux,  will  only  contribute  to  Eq.  (6)  over  area  elements  da 

of  the  control  volume  through  which  an  interface  passes  and  has  a component  of  its  normal 
growth  velocity  vj  along  n.  Identifying  (f)  as  above  and  noting  that  -vi/(u/7/),  we 

see  that 


• n ~ -VI  • n = 


V^M 


VI  n. 


(63) 


where  Eq.  (36)  and  (39)  have  been  used  to  ehminate  e.  For  an  interface  corresponding  to 
Eq.  (35),  one  can  show  from  Eq.  (1)  that  there  is  no  surface  excess  internal  energy  per 
unit  area  relative  to  a sharp  interface  located  at  x = 0.  Therefore,  a = where 

is  the  surface  excess  entropy  per  unit  area  relative  to  a sharp  interface  located  at  x = 0. 
Thus  we  have 


• n da  ~ — vi  • n da. 
V 


(64) 


The  area  over  which  this  term  wiU  contribute  to  Eq.  (6)  is  a thin  ribbon  defined  by  the 
intersection  of  the  diffuse  interface  with  the  area  A of  the  control  volume.  Hence  the  area 
element  da  dir}  csc0,  where  d£  is  an  element  of  length  along  the  ribbon  and  0 is  the 
angle  between  vi  and  n.  This  trigonometric  factor  arises  because  the  diffuse  boundary  is 
not  necessarily  perpendicular  to  the  surface  of  the  control  volume.  Thus  Anally 


J • n da  ~ ^ cot  0 dl,  (65) 

which  is  in  complete  agreement  with  the  analysis  of  Gurtin  for  a sharp  interface  model 
[27,28].  Eq.  (65)  represents  the  rate  at  which  entropy  exits  the  control  volume  by  means  of 
the  passage  of  portions  of  the  interface  out  of  the  surface  of  the  control  volume.  It  makes 
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no  contribution  in  regions  of  that  surface  at  which  the  interface  is  moving  in  a direction 
perpendicular  to  the  normal  of  the  control  volume. 

Penrose  and  Fife  were  able  to  obtain  hanger’s  model  from  a formahsm  that  guarantees 
positive  entropy  production  by  making  the  following  choices  (in  our  notation,  equivalent 
to  their  Eqn.  3.11); 

p(^)  = <t>,  L{T)  = io,  es{T)  = /I  - B/T, 

where  A and  B are  constants.  This  choice  of  p{4>),  in  view  of  Eq.  (27),  would  require  G{^) 
to  depend  on  T,  a contradiction.  Therefore,  for  their  model,  the  values  of  (j)  corresponding 
to  the  bulk  phases  will  depend  on  temperature. 

The  methodology  developed  in  the  present  paper  can  be  extended  to  develop  a phase- 
field  model  for  the  non-isothermal  sohdification  of  an  alloy.  For  a binary  alloy,  for  example, 
it  would  be  necessary  to  specify  the  internal  energy  as  a function  of  T,  (f),  and  the  composi- 
tion c.  This  would  be  a logical  extension  of  the  phase-field  model  for  isothermal  binary  al- 
loy solidification,  developed  for  a lens-type  phase  diagram  in  Ref.  [29].  Another  extension, 
which  will  be  the  subject  of  a forthcoming  paper,  is  to  incorporate  arbitrary  anisotropy  of 
surface  tension  and  interface  kinetics  in  a thermodynamically-consistent  manner. 
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Figure  Captions 


Figure  1.  Sketch  of  the  internal  energies  e£,(r)  and  es{T)  as  functions  of  temperature.  The 
sohd  curves  correspond  to  the  stable  bulk  phases,  whereas  the  dashed  curves  correspond 
to  either  supercooled  liquid  or  superheated  soHd.  The  difference  between  the  curves  at 
any  temperature  is  the  function  L{T),  which  is  equal  to  the  latent  heat  Lq  = L[Tm)  at 
the  melting  point  Tm- 

Figure  2.  Plot  of  the  quantity  —L'{u)  versus  u that  illustrates  (see  Eq.  (60))  the  temperature- 
dependence  of  the  specific  heat  of  the  sohd,  C5,  that  is  necessary  to  obtain  a thermodynamically- 
consistent  equation  for  temperature  corresponding  to  Kobayashi’s  model.  The  plot  is  for 
u = 10,  used  by  Kobayashi  in  his  computations,  and  {Tm  — To)fTM  = 0.2.  For  values  of 
the  undercooling  with  5 < 1,  there  will  be  a rather  large  variation  of  cs  with  temperature, 
even  for  temperatures  near  the  melting  point. 
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pecific  forms  of  a phase-field  function  are  chosen  to  produce  two  models  that  bear  strong 
resemblances  to  the  models  proposed  by  Danger  and  Kobayashi.  Our  models  contain  additional 
nonlinear  functions  of  the  phase  field  that  are  necessary  to  guarantee  thermodynamic 
consistency. 
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